Groundwater.f90 Source File
Simulate quasi 3D groundwater flux and river-groundwater interaction
Source Code
!! Simulate quasi 3D groundwater flux and river-groundwater interaction !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.1 - 12th February 2023 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 14/Sep/2011 | Original code | ! | 1.0 | 22/Feb/2023 | Code rewritten | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description !Simulate quasi 3D groundwater flux and river-groundwater interaction ! !### References ! ! Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic cellular ! automata for groundwater modelling: a first approach. Environmental ! modelling & software, 26(5), 634–643. ! ! Ravazzani, G., Curti, D., Gattinoni, P., Della Valentina, S., ! Fiorucci, A., & Rosso, R.. (2016). Assessing groundwater ! contribution to streamflow of a large alpine river with heat ! tracer methods and hydrological modelling. River research ! and applications, 32(5), 871-884. ! ! Ravazzani, G., Barbero, S., Salandin, A., Senatore, A., ! & Mancini, M.. (2015). An integrated hydrological model ! for assessing climate change impacts on water resources ! of the upper po river basin. Water resources management, ! 29 (4), 1193-1215. ! MODULE Groundwater ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Type Definitions: short, float USE DomainProperties, ONLY : & !imported variables: mask USE GridLib, ONLY : & !imported Type definitions: grid_real, grid_integer, & !imported routines: NewGrid, GridDestroy, & !Imported parameters: NET_CDF USE GridOperations, ONLY : & !Imported routines: GridByIni, CellArea, & !Imported operators: ASSIGNMENT (=) USE IniLib, ONLY : & !Imported definitions: IniList, & !Imported routines: IniOpen, IniClose, & SectionIsPresent, KeyIsPresent, & SubSectionIsPresent, & IniReadReal, IniReadInt USE LogLib, ONLY: & !Imported routines: Catch USE StringManipulation, ONLY : & !Imported routines ToString USE Chronos, ONLY : & !Imported definitions: DateTime , & !Imported variables: timeString, & !Imported operations: OPERATOR (+), & OPERATOR (==), & ASSIGNMENT (=) USE ObservationalNetworks, ONLY : & !Imported routines: ReadMetadata, CopyNetwork, & WriteMetadata, DestroyNetwork, & AssignDataFromGrid, WriteData, & !Imported defined variable: ObservationalNetwork USE Utilities, ONLY : & !Imported routines: GetUnit USE GridStatistics, ONLY : & !imported routines: GetMean USE TableLib, ONLY: & !imported definitions: Table, & !Imported routines: TableNew, TableGetValue IMPLICIT NONE ! Global (i.e. public) Declarations: ! Global TYPE Definitions: TYPE Aquifer TYPE(grid_real) :: top TYPE(grid_real) :: bottom TYPE(grid_real) :: Ks !hydraulic conductivity (m/s) TYPE(grid_real) :: sy !specific yield TYPE(grid_real) :: KsAquitard !hydraulic conductivity of the aquitard TYPE(grid_integer):: domainBC !set domain and boundary conditions type TYPE(grid_real) :: valueBC ! value to be used for boundary condition TYPE(grid_real) :: head0 !head at time t-1 TYPE(grid_real) :: head1 !head at time t END TYPE Aquifer TYPE GroundwaterBasin INTEGER (KIND = short) :: nAquifers !!number of aquifers in the basin TYPE (Aquifer), POINTER :: aquifer (:) END TYPE GroundwaterBasin !Global routines: PUBLIC :: GroundwaterInit PUBLIC :: GroundwaterUpdate PUBLIC :: GroundwaterPointInit PUBLIC :: GroundwaterOutputInit PUBLIC :: GroundwaterOutput PUBLIC :: GroundwaterRiverUpdate !Global variables: TYPE (GroundwaterBasin) :: basin INTEGER (KIND = short) :: dtGroundwater TYPE (grid_real) :: vadoseDepth TYPE (grid_real) :: groundwaterToRiver !!discharge from groundwater to river (m3/s) TYPE (grid_real) :: riverToGroundwater !!discharge from river to groundwater (m3/s) LOGICAL :: riverGroundwaterInteract !Local (private) parameters: INTEGER (KIND = short), PARAMETER, PRIVATE :: BC_DIRICHLET = 2 ! head boundary condition INTEGER (KIND = short), PARAMETER, PRIVATE :: BC_NEUMANN = 3 ! flux boundary condition INTEGER (KIND = short), PARAMETER, PRIVATE :: ACTIVE_CELL = 1 ! active cell within the groundwater domain !local variables TYPE (grid_real), PRIVATE :: transmissivity !!local transmissivity (m2/s) TYPE (grid_real), PRIVATE :: fluxUpward !!flux from lower aquitard (m/s) TYPE (grid_real), PRIVATE :: fluxDownward !!flux from upper aquitard (m/s) TYPE (grid_real), PRIVATE :: riverDem !!reference digital elevation model to compute !!river water surface elevation (m asl) TYPE (grid_real), PRIVATE :: streambedThickness !! streambed thickness (m) TYPE (grid_real), PRIVATE :: streambedConductivity !! streambed conductivity (m/s) TYPE (grid_integer), PRIVATE :: riverGroundwaterID !!id of river reaches that interact with groundwater TYPE (Table), PRIVATE :: riverGroundwaterParameters TYPE (ObservationalNetwork),PRIVATE :: sites INTEGER (KIND = short), ALLOCATABLE, PRIVATE :: fileUnitPointGW (:) INTEGER (KIND = short), PRIVATE :: fileUnitOutGW TYPE (DateTime), PRIVATE :: timePointExport REAL (KIND = float), PRIVATE :: volumeRecharge !! recharge volume (m3) REAL (KIND = float), PRIVATE :: volumeResidual !! balance error (m3) REAL (KIND = float), PRIVATE :: volumeNeumann !!volume from Neumann boundary condition (m3) REAL (KIND = float), PRIVATE :: volumeDirichlet !!volume from Dirichlet boundary condition (m3) REAL (KIND = float), PRIVATE :: volumeChange !! change of stored water volume in groundwater (m3) REAL (KIND = float), PRIVATE :: volumeRiverToGroundwater !! volume from river to groundwater (m3) REAL (KIND = float), PRIVATE :: volumeGroundwaterToRiver !! volume from groundwater to river (m3) !local (private) routines: PRIVATE :: GroundwaterParameterUpdate PRIVATE :: GroundwaterPointExport PRIVATE :: GroundwaterRiverInit !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize groundwater SUBROUTINE GroundwaterInit & ! (inifile) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information !local declarations: TYPE (IniList) :: iniDB INTEGER (KIND = short) :: k , i, j CHARACTER (LEN = 100) :: string REAL (KIND = float) :: scalar !-----------------------------------end of declarations------------------------ !open and read configuration file CALL IniOpen (inifile, iniDB) !read the number of aquifers IF ( KeyIsPresent ('aquifers', iniDB ) ) THEN basin % naquifers = IniReadInt ( 'aquifers', iniDB ) ELSE CALL Catch ('error', 'Groundwater', & 'aquifers not found in configuration file' ) END IF !allocate aquifers ALLOCATE ( basin % aquifer ( basin % naquifers ) ) !load parameters DO k = 1, basin % naquifers string = 'aquifer_' // ToString (k) IF (SectionIsPresent(TRIM(string), iniDB)) THEN !load domain and boundary condition id IF (SubSectionisPresent (subsection = 'boundary-condition-id', & section = TRIM(string), iniDB = iniDB) ) THEN CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % domainBC, & section = TRIM(string), & subsection = 'boundary-condition-id') ELSE CALL Catch ('error', 'Groundwater', & 'missing boundary-condition-id in configuration file: ', & argument = TRIM(string) ) END IF !load boundary condition value IF (SubSectionisPresent (subsection = 'boundary-condition-value', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'boundary-condition-value' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'boundary-condition-value' ) CALL NewGrid (basin % aquifer (k) % valueBC, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % valueBC, & section = TRIM(string), & subsection = 'boundary-condition-value') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing boundary-condition-value in configuration file: ', & argument = TRIM(string) ) END IF !load initial head IF (SubSectionisPresent (subsection = 'initial-head', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'initial-head' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'initial-head' ) CALL NewGrid (basin % aquifer (k) % head0, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % head0, & section = TRIM(string), & subsection = 'initial-head') END IF CALL NewGrid (basin % aquifer (k) % head1, & basin % aquifer (k) % domainBC, 0.) basin % aquifer (k) % head1 = basin % aquifer (k) % head0 !head boundary condition overlay DO i = 1, basin % aquifer (k) % domainBC % idim DO j = 1, basin % aquifer (k) % domainBC % jdim IF ( basin % aquifer (k) % domainBC % mat (i,j) == & BC_DIRICHLET ) THEN basin % aquifer (k) % head0 % mat (i,j) = & basin % aquifer (k) % valueBC % mat (i,j) basin % aquifer (k) % head1 % mat (i,j) = & basin % aquifer (k) % valueBC % mat (i,j) END IF END DO END DO ELSE CALL Catch ('error', 'Groundwater', & 'missing initial-condition-head in configuration file: ', & argument = TRIM(string) ) END IF !load top elevation IF (SubSectionisPresent (subsection = 'top-elevation', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'top-elevation' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'top-elevation' ) CALL NewGrid (basin % aquifer (k) % top, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % top, & section = TRIM(string), & subsection = 'top-elevation') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing top-elevation in configuration file: ', & argument = TRIM(string) ) END IF !load bottom elevation IF (SubSectionisPresent (subsection = 'bottom-elevation', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'bottom-elevation' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'bottom-elevation' ) CALL NewGrid (basin % aquifer (k) % bottom, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % bottom, & section = TRIM(string), & subsection = 'bottom-elevation') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing bottom-elevation in configuration file: ', & argument = TRIM(string) ) END IF !load hydraulic conductivity IF (SubSectionisPresent (subsection = 'hydraulic-conductivity', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'hydraulic-conductivity' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'hydraulic-conductivity' ) CALL NewGrid (basin % aquifer (k) % Ks, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % Ks, & section = TRIM(string), & subsection = 'hydraulic-conductivity') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing hydraulic-conductivity in configuration file: ', & argument = TRIM(string) ) END IF !load specific yield IF (SubSectionisPresent (subsection = 'specific-yield', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'specific-yield' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'specific-yield' ) CALL NewGrid (basin % aquifer (k) % sy, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % sy, & section = TRIM(string), & subsection = 'specific-yield') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing specific-yield in configuration file: ', & argument = TRIM(string) ) END IF ! aquitard hydraulic conductivity IF ( k < basin % naquifers ) THEN IF (SubSectionisPresent (subsection = 'aquitard-conductivity', & section = TRIM(string), iniDB = iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, TRIM(string), & 'aquitard-conductivity' ) ) THEN scalar = IniReadReal ('scalar', iniDB, & TRIM(string), 'aquitard-conductivity' ) CALL NewGrid (basin % aquifer (k) % KsAquitard, & basin % aquifer (k) % domainBC, scalar) ELSE CALL GridByIni ( ini = iniDB, & grid = basin % aquifer (k) % KsAquitard, & section = TRIM(string), & subsection = 'aquitard-conductivity') END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing aquitard-conductivity in configuration file: ', & argument = TRIM(string) ) END IF END IF ELSE CALL Catch ('error', 'Groundwater', & 'missing section in configuration file: ', & argument = TRIM(string) ) END IF END DO !allocate variables CALL NewGrid ( transmissivity, basin % aquifer (1) % domainBC, 0. ) CALL NewGrid ( fluxUpward, basin % aquifer (1) % domainBC, 0. ) CALL NewGrid ( fluxDownward, basin % aquifer (1) % domainBC, 0. ) CALL NewGrid ( vadoseDepth, basin % aquifer (1) % domainBC, 0. ) !check river-groundwater interaction IF ( SectionIsPresent ( 'river-groundwater', iniDB) ) THEN riverGroundwaterInteract = .TRUE. CALL GroundwaterRiverInit ( inifile ) ELSE riverGroundwaterInteract = .FALSE. END IF !destroy iniDB CALL IniClose (IniDB) END SUBROUTINE GroundwaterInit !============================================================================== !| Description: ! update groundwater head with Darcy law and mass conservation ! in two dimensions, with a macrocopic cellular automata approach !``` ! +-----+ ! | N | ! +-----+-----+-----+ ! | W | C | E | ! +-----+-----+-----+ ! | S | ! +-----+ !``` ! References: ! Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic cellular ! automata for groundwater modelling: a first approach. Environmental ! modelling & software, 26(5), 634–643. ! SUBROUTINE GroundwaterUpdate & ! ( time, hillslopeFlux, percolation, capflux ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time TYPE (grid_real), INTENT(IN) :: hillslopeFlux !!flux from hillslope (m3/s) TYPE (grid_real), INTENT(IN) :: percolation !! depp infiltration from soil balance (m/s) TYPE (grid_real), INTENT(IN) :: capflux !! capillary flux from soil balance (m/s) !local declarations: INTEGER (KIND = short) :: i, j, k REAL (KIND = float) :: rechargeQ !!vertical recharge (m3/s) REAL (KIND = float) :: transNC !!mean transmissivity North-Centre (m2/s) REAL (KIND = float) :: transEC !!mean transmissivity East-Centre (m2/s) REAL (KIND = float) :: transSC !!mean transmissivity South-Centre (m2/s) REAL (KIND = float) :: transWC !!mean transmissivity West-Centre (m2/s) REAL (KIND = float) :: fluxNC !!flux North-Centre (m3/s) REAL (KIND = float) :: fluxEC !!flux East-Centre (m3/s) REAL (KIND = float) :: fluxSC !!flux South-Centre (m3/s) REAL (KIND = float) :: fluxWC !!flux West-Centre (m3/s) REAL (KIND = float) :: fluxNESW !! total flux from the 4 directions (m3/s) REAL (KIND = float) :: areaOfCell !! cell area (m2) REAL (KIND = float) :: rToG !!river to groundwater discharge (m3/s) REAL (KIND = float) :: gToR !!groundwater to river discharge (m3/s) !-------------------------------------end of declarations---------------------- !reset counter variables volumeNeumann = 0. volumeDirichlet = 0. volumeRecharge = 0. volumeChange = 0. volumeResidual = 0. !update boundary condition CALL GroundwaterParameterUpdate (time) !save head1 into head0 DO k = 1, basin % nAquifers basin % aquifer (k) % head0 = basin % aquifer (k) % head1 END DO !update flux from hillslope of first (surface) aquifer DO i = 1, basin % aquifer (1) % domainBC % idim DO j = 1, basin % aquifer (1) % domainBC % jdim IF ( basin % aquifer (1) % domainBC % mat (i,j) == & BC_NEUMANN ) THEN basin % aquifer (1) % valueBC % mat (i,j) = & hillslopeFlux % mat (i,j) END IF END DO END DO !update aquifers head DO k = 1, basin % nAquifers !compute local transmissivity DO i = 1, basin % aquifer (k) % domainBC % idim DO j = 1, basin % aquifer (k) % domainBC % jdim IF ( basin % aquifer (k) % domainBC % mat(i,j) /= & basin % aquifer (k) % domainBC % nodata ) THEN !transmissivity: for confined aquifer = ks * aquifer depth transmissivity % mat (i,j) = MIN ( & basin % aquifer (k) % Ks % mat (i,j) * & ( basin % aquifer (k) % head0 % mat (i,j) - & basin % aquifer (k) % bottom % mat (i,j) ) , & basin % aquifer (k) % Ks % mat (i,j) * & ( basin % aquifer (k) % top % mat (i,j) - & basin % aquifer (k) % bottom % mat (i,j) ) & ) IF ( transmissivity % mat (i,j) <= 0. ) THEN transmissivity % mat (i,j) = 0.01 END IF END IF END DO END DO !update head1 DO i = 1, basin % aquifer (k) % domainBC % idim DO j = 1, basin % aquifer (k) % domainBC % jdim IF ( basin % aquifer (k) % domainBC % mat(i,j) == & ACTIVE_CELL ) THEN !harmonic mean transmissivity transNC = 2.0 * transmissivity % mat(i-1,j) * & transmissivity % mat(i,j) / & ( transmissivity %mat(i-1,j) + & transmissivity %mat(i,j) ) transSC = 2.0 * transmissivity % mat(i+1,j) * & transmissivity % mat(i,j) / & ( transmissivity % mat(i+1,j) + & transmissivity % mat(i,j) ) transEC = 2.0 * transmissivity % mat(i,j+1) * & transmissivity % mat(i,j) / & ( transmissivity % mat(i,j+1) + & transmissivity % mat(i,j) ) transWC = 2.0 * transmissivity % mat(i,j-1) * & transmissivity % mat(i,j) / & ( transmissivity % mat(i,j-1) + & transmissivity % mat(i,j) ) !flux from North IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == & BC_NEUMANN ) THEN fluxNC = basin % aquifer (k) % valueBC % mat(i-1,j) volumeNeumann = volumeNeumann + fluxNC * dtGroundwater ELSE fluxNC = transNC * & ( basin % aquifer (k) % head0 % mat (i-1,j) - & basin % aquifer (k) % head0 % mat (i,j) ) IF ( basin % aquifer (k) % domainBC % mat(i-1,j) == & BC_DIRICHLET ) THEN volumeDirichlet = volumeDirichlet + & fluxNC * dtGroundwater END IF END IF !flux from East IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == & BC_NEUMANN ) THEN fluxEC = basin % aquifer (k) % valueBC % mat(i,j+1) volumeNeumann = volumeNeumann + fluxEC * dtGroundwater ELSE fluxEC = transEC * & ( basin % aquifer (k) % head0 % mat (i,j+1) - & basin % aquifer (k) % head0 % mat (i,j) ) IF ( basin % aquifer (k) % domainBC % mat(i,j+1) == & BC_DIRICHLET ) THEN volumeDirichlet = volumeDirichlet + & fluxEC * dtGroundwater END IF END IF !flux from South IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == & BC_NEUMANN ) THEN fluxSC = basin % aquifer (k) % valueBC % mat(i+1,j) volumeNeumann = volumeNeumann + fluxSC * dtGroundwater ELSE fluxSC = transSC * & ( basin % aquifer (k) % head0 % mat (i+1,j) - & basin % aquifer (k) % head0 % mat (i,j) ) IF ( basin % aquifer (k) % domainBC % mat(i+1,j) == & BC_DIRICHLET ) THEN volumeDirichlet = volumeDirichlet + & fluxSC * dtGroundwater END IF END IF !flux from West IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == & BC_NEUMANN ) THEN fluxWC = basin % aquifer (k) % valueBC % mat(i,j-1) volumeNeumann = volumeNeumann + fluxWC * dtGroundwater ELSE fluxWC = transWC * & ( basin % aquifer (k) % head0 % mat (i,j-1) - & basin % aquifer (k) % head0 % mat (i,j) ) IF ( basin % aquifer (k) % domainBC % mat(i,j-1) == & BC_DIRICHLET ) THEN volumeDirichlet = volumeDirichlet + & fluxWC * dtGroundwater END IF END IF !total flux fluxNESW = fluxNC + fluxEC + fluxSC + fluxWC !area of the cell areaOfCell = CellArea (transmissivity, i, j) !vertical recharge IF ( ALLOCATED (percolation % mat ) ) THEN IF ( k == 1 .AND. percolation % mat (i,j) /= & percolation % nodata ) THEN rechargeQ = ( percolation % mat (i,j) - & capflux % mat (i,j) ) * areaOfCell volumeRecharge = volumeRecharge + & rechargeQ * dtGroundwater ELSE rechargeQ = 0. END IF ELSE rechargeQ = 0. END IF !upward flux from lower aquifer IF (k /= basin % nAquifers ) THEN fluxUpward % mat(i,j) = & basin % aquifer (k) % ksAquitard % mat(i,j) * & ( ( basin % aquifer (k+1) % head0 % mat(i,j) - & basin % aquifer (k) % head0 % mat(i,j) ) / & ( basin % aquifer (k) % bottom % mat(i,j) - & basin % aquifer (k+1) % top % mat(i,j) ) ) ELSE fluxUpward % mat (i,j) = 0.0 END IF !update head1 IF ( k == 1) THEN fluxDownward % mat (i,j) = 0. IF ( riverGroundwaterInteract ) THEN rToG = riverToGroundwater % mat (i,j) gToR = - groundwaterToRiver % mat (i,j) ELSE rToG = 0. gToR = 0. END IF ELSE rToG = 0. gToR = 0. END IF basin % aquifer (k) % head1 % mat(i,j) = & basin % aquifer (k) % head0 % mat(i,j) + & ( fluxNESW + rechargeQ + rToG + gToR ) * & dtGroundwater / areaOfCell / & basin % aquifer (k) % sy % mat(i,j) + & ( fluxUpward % mat(i,j) + fluxDownward % mat(i,j) ) * & dtGroundwater / basin % aquifer (k) % sy % mat(i,j) IF ( ISNAN ( basin % aquifer (k) % head1 % mat(i,j) ) ) THEN basin % aquifer (k) % head1 % mat(i,j) = & basin % aquifer (k) % head0 % mat(i,j) volumeChange = 0. ELSE !volume change volumeChange = volumeChange + & ( basin % aquifer (k) % head1 % mat(i,j) - & basin % aquifer (k) % head0 % mat(i,j) ) * & basin % aquifer (k) % sy % mat(i,j) * areaOfCell END IF !swap fluxUpward and fluxDownward fluxDownward % mat (i,j) = - fluxUpward % mat (i,j) !update vadose zone depth IF ( k == 1) THEN vadoseDepth % mat (i,j) = & basin % aquifer (k) % top % mat(i,j) - & basin % aquifer (k) % head1 % mat(i,j) END IF END IF END DO END DO END DO !compute balance error volumeResidual = volumeRecharge + volumeDirichlet + & volumeNeumann + volumeRiverToGroundwater - & volumeGroundwaterToRiver - volumeChange !write point data IF ( time == timePointExport ) THEN CALL GroundwaterPointExport ( time ) END IF RETURN END SUBROUTINE GroundwaterUpdate !============================================================================== !| Description: ! update boundary condition map that change in time SUBROUTINE GroundwaterParameterUpdate & ! (time) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: CHARACTER (LEN = 300) :: filename CHARACTER (LEN = 300) :: varname INTEGER (KIND = short) :: k, i, j !------------------------------end of declarations----------------------------- !boundary condition DO k = 1, basin % naquifers IF ( time == basin % aquifer (k) % valueBC % next_time ) THEN !destroy current grid filename = basin % aquifer (k) % valueBC % file_name varname = basin % aquifer (k) % valueBC % var_name CALL GridDestroy (basin % aquifer (k) % valueBC ) !read new grid in netcdf file CALL NewGrid (basin % aquifer (k) % valueBC, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) !head boundary condition overlay DO i = 1, basin % aquifer (k) % domainBC % idim DO j = 1, basin % aquifer (k) % domainBC % jdim IF ( basin % aquifer (k) % domainBC % mat (i,j) == & BC_DIRICHLET ) THEN basin % aquifer (k) % head0 % mat (i,j) = & basin % aquifer (k) % valueBC % mat (i,j) basin % aquifer (k) % head1 % mat (i,j) = & basin % aquifer (k) % valueBC % mat (i,j) END IF END DO END DO END IF END DO RETURN END SUBROUTINE GroundwaterParameterUpdate !============================================================================== !| Description: ! Initialize export of point site data SUBROUTINE GroundwaterPointInit & ! ( pointfile, path_out, time ) IMPLICIT NONE !Arguments with intent (in): CHARACTER (LEN = *), INTENT(IN) :: pointfile !!file containing coordinate of points CHARACTER (LEN = *), INTENT(IN) :: path_out !!path of output folder TYPE (DateTime), INTENT(IN) :: time !!start time !local declarations INTEGER (KIND = short) :: err_io INTEGER (KIND = short) :: fileunit INTEGER (KIND = short) :: k !-------------------------end of declarations---------------------------------- timePointExport = time fileunit = GetUnit () OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), & status='OLD', iostat = err_io ) IF ( err_io /= 0 ) THEN !file does not exist CALL Catch ('error', 'Groundwater', 'out point file does not exist') END IF !Read metadata CALL ReadMetadata (fileunit, sites) !check dt IF (.NOT. ( MOD ( sites % timeIncrement, dtGroundwater ) == 0 ) ) THEN CALL Catch ('error', 'Groundwater', 'dt out sites must be multiple of dtGroundwater ') END IF CLOSE ( fileunit ) sites % description = 'groundwater head data exported from FEST' sites % unit = 'm' sites % offsetZ = 0. !allocate file unit ALLOCATE ( fileUnitPointGW ( basin % nAquifers ) ) !open and initialize files DO k = 1, basin % nAquifers fileUnitPointGW (k) = GetUnit () OPEN ( unit = fileUnitPointGW (k), & file = TRIM(path_out) // 'point_aquifer_' // TRIM(ToString (k)) // '.fts' ) CALL WriteMetadata ( network = sites, fileunit = fileUnitPointGW (k) ) CALL WriteData (sites, fileUnitPointGW (k), .TRUE.) END DO RETURN END SUBROUTINE GroundwaterPointInit !============================================================================== !| Description: ! Export of point site data SUBROUTINE GroundwaterPointExport & ! ( time ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: INTEGER (KIND = short) :: k !-------------------------end of declarations---------------------------------- !set current time sites % time = time DO k = 1, basin % nAquifers !populate data CALL AssignDataFromGrid ( basin % aquifer (k) % head1, sites ) !write data CALL WriteData ( sites, fileUnitPointGW (k) ) END DO timePointExport = timePointExport + sites % timeIncrement RETURN END SUBROUTINE GroundwaterPointExport !============================================================================== !| Description: ! Initialize files for exporting output results SUBROUTINE GroundwaterOutputInit & ! ( path_out ) IMPLICIT NONE !Arguments with intent (in): CHARACTER (LEN = *), INTENT(IN) :: path_out !!path of output folder !local declarations REAL (KIND = float) :: area !! (km2) INTEGER (KIND = short) :: i, j INTEGER (KIND = short) :: nvars !----------------------------------end of declarations------------------------- !compute groundwater extent area area = 0. DO i = 1, basin % aquifer (1) % domainBC % idim DO j = 1, basin % aquifer (1) % domainBC % jdim IF (basin % aquifer (1) % domainBC % mat (i,j) == ACTIVE_CELL ) THEN area = area + CellArea ( basin % aquifer (1) % domainBC, i, j ) END IF END DO END DO area = area / 1000000. !conversion m2 -> km2 !number of variables to write in output file nvars = basin % nAquifers + 5 IF ( riverGroundwaterInteract ) THEN nvars = nvars + 2 END IF !open and initialize file fileUnitOutGW = GetUnit () OPEN ( unit = fileUnitOutGW, file = TRIM(path_out) // 'aquifer.out' ) WRITE (fileUnitOutGW,'(a)') 'Output variables of groundwater simulation' WRITE (fileUnitOutGW,fmt = '(a, F12.2)') 'extent area (km2): ', area WRITE (fileUnitOutGW,fmt = '(a, I3)') 'number of variables: ', nvars WRITE (fileUnitOutGW,*) WRITE (fileUnitOutGW,*) WRITE (fileUnitOutGW,'(a)') 'data' WRITE (fileUnitOutGW,'(a)', advance='no') 'DateTime ' DO i = 1, basin % nAquifers WRITE (fileUnitOutGW,'(a)', advance='no') 'head_aquifer_' // & TRIM(ToString(i)) // '_m' END DO IF ( riverGroundwaterInteract ) THEN WRITE (fileUnitOutGW,'(a)') ' recharge_m3 volume_neumann_m3 & volume_dirichlet_m3 volume_change_m3 residual_m3 & river_to_groundwater_m3 groundwater_to_river_m3 ' ELSE WRITE (fileUnitOutGW,'(a)') ' recharge_m3 volume_neumann_m3 & volume_dirichlet_m3 volume_change_m3 residual_m3 ' END IF RETURN END SUBROUTINE GroundwaterOutputInit !============================================================================== !! Description: !! Write results on output file SUBROUTINE GroundwaterOutput & ! ( time ) IMPLICIT NONE !Arguments with intent (in): TYPE (DateTime), INTENT(IN) :: time !! current time !local variables: INTEGER (KIND = short) :: k REAL (KIND = float) :: meanHead !----------------------------end of declarations------------------------------- timeString = time WRITE (fileUnitOutGW,'(a)',advance='no') timeString DO k = 1, basin % nAquifers meanHead = GetMean ( basin % aquifer (k) % head1, & maskInteger = basin % aquifer (k) % domainBC ) WRITE (fileUnitOutGW,'(6x,E12.5)',advance='no') meanHead END DO IF ( riverGroundwaterInteract ) THEN WRITE (fileUnitOutGW,'(7(5x,E12.5))') volumeRecharge, volumeNeumann, & volumeDirichlet, volumeChange, volumeResidual, & volumeRiverToGroundwater, volumeGroundwaterToRiver ELSE WRITE (fileUnitOutGW,'(5(5x,E12.5))') volumeRecharge, volumeNeumann, & volumeDirichlet, volumeChange, volumeResidual END IF RETURN END SUBROUTINE GroundwaterOutput !============================================================================== !| Description: ! Configure river-groundwater interaction SUBROUTINE GroundwaterRiverInit & ! ( inifile ) IMPLICIT NONE !Arguments with intent (in): CHARACTER (LEN = *), INTENT(IN) :: inifile !! configuration file !local variables: TYPE (IniList) :: iniDB INTEGER (KIND = short) :: i,j !----------------------------end of declarations------------------------------- !open and read configuration file CALL IniOpen ( inifile, iniDB ) !load river id IF (SubSectionisPresent (subsection = 'river-id', & section = 'river-groundwater', iniDB = iniDB) ) THEN CALL GridByIni ( ini = iniDB, & grid = riverGroundwaterID, & section = 'river-groundwater', & subsection = 'river-id') ELSE CALL Catch ('error', 'Groundwater', & 'missing river-id in configuration file' ) END IF !load river dem IF (SubSectionisPresent (subsection = 'river-dem', & section = 'river-groundwater', iniDB = iniDB) ) THEN CALL GridByIni ( ini = iniDB, & grid = riverDem, & section = 'river-groundwater', & subsection = 'river-dem') ELSE CALL Catch ('error', 'Groundwater', & 'missing river-dem in configuration file' ) END IF !exchange parameters CALL TableNew ( file = inifile, tab = riverGroundwaterParameters, & id = 'river-groundwater') !allocate exchange flux maps CALL NewGrid ( riverToGroundwater, riverGroundwaterID, 0. ) CALL NewGrid ( groundwaterToRiver, riverGroundwaterID, 0. ) !allocate streambed parameter maps CALL NewGrid ( streambedThickness, riverGroundwaterID, 0. ) CALL NewGrid ( streambedConductivity, riverGroundwaterID, 0. ) !populate streambed parameter maps DO i = 1, riverGroundwaterID % idim DO j = 1, riverGroundwaterID % jdim IF ( riverGroundwaterID % mat (i,j) /= & riverGroundwaterID % nodata ) THEN CALL TableGetValue ( & valueIn = REAL(riverGroundwaterID % mat (i,j)),& tab = riverGroundwaterParameters, & keyIn = 'id', keyOut ='streambed-conductivity', & match = 'exact', valueOut = streambedConductivity % mat (i,j) ) CALL TableGetValue ( & valueIn = REAL(riverGroundwaterID % mat (i,j)),& tab = riverGroundwaterParameters, & keyIn = 'id', keyOut ='streambed-thickness', & match = 'exact', valueOut = streambedThickness % mat (i,j) ) END IF END DO END DO !free db CALL IniClose ( iniDB ) RETURN END SUBROUTINE GroundwaterRiverInit !============================================================================== !| Description: ! Update river-groundwater exchange fluxes SUBROUTINE GroundwaterRiverUpdate & ! ( waterDepth, topWidth ) IMPLICIT NONE !Arguments with intent (in): TYPE (grid_real), INTENT(IN) :: waterDepth !!river water depth (m) TYPE (grid_real), INTENT(IN) :: topWidth !!river top width (m) !local variables: INTEGER (KIND = short) :: i, j REAL (KIND = float) :: riverWSE !river water surface elevation (m asl) !TYPE (grid_real) :: waterTable !!aquifer water table (m) REAL (KIND = float) :: waterTable REAL (KIND = float) :: area !-----------------------------end of declarations------------------------------ volumeRiverToGroundwater = 0. volumeGroundwaterToRiver = 0. !waterTable = basin % aquifer (1) % head1 DO i = 1, riverGroundwaterID % idim DO j = 1, riverGroundwaterID % jdim IF ( riverGroundwaterID % mat (i,j) /= & riverGroundwaterID % nodata ) THEN waterTable = basin % aquifer (1) % head1 % mat (i,j) riverWSE = waterDepth % mat (i,j) + riverDem % mat (i,j) area = CellArea (riverGroundwaterID, i, j) IF ( waterTable > riverWSE ) THEN groundwaterToRiver % mat (i,j) = & ( waterTable - riverWSE ) * & !head difference streambedConductivity % mat (i,j) / & !conductivity streambedThickness % mat (i,j) * & !thickness topWidth % mat (i,j) * & !river width area ** 0.5 !cell size riverToGroundwater % mat (i,j) = 0. ELSE IF ( waterTable < riverWSE .AND. & waterTable > riverDem % mat (i,j) ) THEN ! when waterdepth > 10 cm, compute river discharge toward groundwater IF ( waterDepth % mat (i,j) > 0.10 ) THEN riverToGroundwater % mat(i,j) = & ( riverWSE - waterTable ) * & !head difference streambedConductivity % mat (i,j) / & !conductivity streambedThickness % mat (i,j) * & !thickness topWidth % mat (i,j) * & !river width area ** 0.5 !cell size groundwaterToRiver % mat (i,j) = 0. ELSE riverToGroundwater % mat(i,j) = 0. END IF ELSE IF ( waterTable < riverDem % mat (i,j) ) THEN ! when waterdepth > 10 cm, compute river discharge toward groundwater IF ( waterDepth % mat (i,j) > 0.10 ) THEN riverToGroundwater % mat(i,j) = & waterDepth % mat (i,j) * & !head streambedConductivity % mat (i,j) / & !conductivity streambedThickness % mat (i,j) * & !thickness topWidth % mat (i,j) * & !river width area ** 0.5 !cell size groundwaterToRiver % mat (i,j) = 0. ELSE riverToGroundwater % mat(i,j) = 0. END IF END IF volumeRiverToGroundwater = volumeRiverToGroundwater + & riverToGroundwater % mat (i,j) * & dtGroundwater volumeGroundwaterToRiver = volumeGroundwaterToRiver + & groundwaterToRiver % mat (i,j) * & dtGroundwater END IF END DO END DO RETURN END SUBROUTINE GroundwaterRiverUpdate END MODULE Groundwater